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ABSTRACT. We consider a model for an epidemic in a population that occupies geographically 
distinct locations. The disease is spread within subpopulations by contacts between infective 
and susceptible individuals, and is spread between subpopulations by the migration of infected 
individuals. We show how susceptible individuals can act collectively to limit the spread of 
disease during the initial phase of an epidemic, by specifying the distribution that minimises 
the growth rate of the epidemic when the infectives are migrating so as to maximise the growth 
rate. We also give an explicit strategy that minimises the basic reproduction number, which is 
also shown be optimal in terms of the probability of extinction and total size of the epidemic. 
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1. Introduction 


Recently, a number of papers have been devoted to the issue of controlling disease outbreaks. 
Typical mechanisms for control involve treatments which speed recovery [25|, |28j, culling of 

reducing the density of disease vectors [23[], vaccination programs 


infected individuals 

[hi 2o| and quarantine [28j. When the population has some spatial structure, migration also 
plays an important role in disease spread and provides a further control mechanism. 

A common approach to incorporating spatial structure in epidemic modelling is to impose 
a metapopulation structure on the population [see t], 13, 14, 17, for example]. In a metapop¬ 
ulation, the population is divided into a number of subpopulations occupying geographically 
distinct locations. The disease is spread within a subpopulation by contacts between infective 
and susceptible individuals and is spread between subpopulations by the migration of infected 
individuals. 

The effect of migration rates on disease spread in metapopulations has been investigated in a 

number of papers. Due to the complexity of these models, control strategies are often based on 
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minimising the basic reproduction number Rq. Studying a multi-patch frequency dependent SIS 


model, Allen et al [2] note that the rapid movement of infective individuals can lead to disease 
extinction in low risk environments. Furthermore, they conjecture that Rq is a decreasing 


function of the diffusion rate for infective individuals. Hsieh et al |18l . Theorem 4.2] note a 
similar result for their two-patch SEIRP model and a similar phenomena has been observed 
in population models with spatially heterogeneous environments [10]. However, Gao and R.uan 


12l . section 4] have shown that for other models the dependence of Rq on migration rates can 


be more complex. To investigate the effect of the migration rates on other quantities sue! 
the number of infected individuals, numerical methods are generally required [for example 


i as 


29]. 


In this paper, we examine how susceptible individuals can act collectively to limit the spread 
of disease during the initial phase of an epidemic. More specifically, we consider how susceptible 
individuals can distribute themselves in the metapopulation in a way that minimises the growth 
of the epidemic when the infectives migrate so as to maximise the growth. By formulating the 
problem as a minimax optimisation and focusing on the susceptible individuals, we avoid the 
need to distinguish between infected and susceptible individuals when applying controls to the 
population. This is advantageous as identification of infected individuals can be problematic due 
to factors such as delays in the onset of symptoms, asymptomatic carriers and costs associated 
with testing. Furthermore, acute disease can have a significant effect on the behaviour of animals 


15[. This is particularly true for certain parasitic diseases where the parasite attempts to force 


0 . 


the host to act in a manner which assists the propagation of the parasite 

In Section 2 we give our main results. Instead of using an ODE model for the epidemic as was 
done in the papers cited above, our analysis is based on a branching process model. Branching 
processes are known to provide a good approximation to the standard SIR and SIS Markov 


chain models when the number of infectives is initially small Q]. Using this model, we are able 
to give an explicit strategy that minimises the expected rate of growth under a certain condition 
on the recovery and infection rates. We also give an explicit strategy that minimises the basic 
reproduction number which does not require this extra condition. This later strategy is shown 
to also be optimal in terms of the probability of extinction and total size of the epidemic. In 
Section 3, the problem of minimising the expected growth rate is investigated numerically. The 
paper concludes with a discussion of how the results depend on contact rates and how they 
relate to ODE models. 
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2. Minimising disease spread in the initial stages 

Consider a closed population of size N divided into m groups such that at time t group i 
contains Xi(t) susceptibles and ft(t) infectives. Each individual, conditional on its disease 
status, moves independently between groups according to an irreducible Markov process on 
{1,..., m} with transition rate matrix R if it is susceptible and transition rate matrix Q if it is 
infected. The epidemic evolves as a Markov process. Contacts between individuals in the same 
group are assumed to be density dependent jj3]. More precisely, a pair of individuals in group 
i makes contact at the points of a Poisson process of rate /3i/N with contacts between distinct 
pairs of individuals being mutually independent. It is assumed that contact between an infective 
and a susceptible results in the infection of the susceptible. An infected individual in group i 
recovers with immunity at a rate 7 j. Since we are primarily concerned with the initial phase 
of the epidemic, our conclusions remain valid for epidemics where individuals recover without 
immunity. 

In the absence of infective individuals, the entirely susceptible population evolves following 
a closed (linear) migration process with per-capita migration rates R. If the population is in 
equilibrium, then the probability that an individual is in group i is given by 77 where 7 r is the 
unique solution to nR = 0 subject to the constraint 7rl = 1. 

We consider the spread of the disease from a small number of initial infective individuals. 
Clancy [ 7 , Theorem 2.1] shows that, when N is large, the epidemic can be approximated by 
a multi-type branching process. Assuming the susceptible population is in equilibrium, the 
branching process for the number of infective individuals is given by 


(n,. 


->(■• 

• > ft + 1,.. 

■,Yj- 1,...) 

at rate QjiYj, 

(2.1) 

(n,. 

■ ■,Y m ) 

->(■■ 

■ ,Yi + 1, .. 

•) 

at rate ft77ft, 

(2.2) 

m,- 

■ ■,Y m ) 


-,ft-l,-- 

•) 

at rate 7 * ft- 

(2.3) 


Note that the branching process depends on R only through the equilibrium distribution 7r. 

Suppose that the susceptible population aims to minimise some quantity /( ir, Q), calculated 
from the branching process determined by Let S denote the relative interior of the 

(m — l)-simplex and let Q be the set of irreducible migration rate matrices. Without imposing 
any constraints on the movements of the infectives, the susceptible population can choose 7r 
such that, for any e > 0, a value no larger than inftgs sup^gQ f(ir, Q) + e is attained. On the 
other hand, the infectives can migrate in such a way that, for any e > 0, a value no smaller 
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than supg g Q inf 7re< s/(7r, Q) — e is attained. In general, 


sup inf /(7r, Q) < inf sup f(n,Q) 

QeQir&S QeQ 


32l . Lemma in section 1.2.2], A pair (tt*,Q*) £ S x Q such that 




for all 7r £ S and all Q £ Q is called a saddle point for /. If a saddle point exists, then 


min sup /(7r, Q) = max inf /(tt, Q ) 
ireS q £ q QeQneS 


32, Theorem in section 1.3.4]. The susceptibles can attain this value by distributing themselves 
amongst the groups according to n*. When a saddle point for / does not exist, there may still 
be an e-saddle point, that is, for every e > 0 there exists a pair ( tt € , Q e ) £ S x Q such that 


/O', Q)-e< /(tt £ , Q £ ) < /(tt, Q e ) + e, 


for all 7r € S and all Q £ Q. The existence of an e-saddle point implies that 


inf sup /(7 t,Q) = sup inf /(7 t,Q) = lim/(7r e , Q e ) 


32, Theorem in section 2.2.5]. In the following, we determine the (e)-saddle points for four 


quantities derived from the branching process dziii-ra- 

As mentioned in the introduction, this formulation avoids the need to distinguish between 
susceptible and infected individuals in the application of controls. To illustrate this point, 
suppose that susceptible individuals normally move between groups following a Markov process 
with migration rate matrix R. The optimal distribution of susceptibles n* can be obtained by 
border controls where a migrating individual from group j going to group i is given admittance 
with probability pji and otherwise returned to group j. Detailed balance equations show that 
the optimal distribution for susceptibles is obtained if the admittance probabilities satisfy 


for all i,j. Although the border controls will have an effect on the migration rate of in¬ 
fected individuals if they are applied to the population as a whole, the optimal distribution 
for susceptible individuals ensures that the growth of the epidemic can be no greater than 
min^s sup QeQ /(tt, Q). 
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2.1. Minimising the expected growth rate. Let {M(t);t > 0} be the mean matrix semi¬ 
group 

Mij(t ) = E(Y)(t) | y r (0) = S ri , r = 1,2,... ,m), 

for i,j G {1, 2,... , m} and where the Sij are Kronecker deltas. For a vector a G M m , let 
diag(a) denote the m x m diagonal matrix with diag {a)a = Oj, i = 1,..., m. The mean matrix 
semigroup has the infinitesimal generator 


yl(7r, Q) = diag(/3)diag(?r) - diag( 7 ) + Q 

and M(t ) = exp(A(n,Q)t) for t > 0 [fi]. By Seneta 30, Theorem 2.7], if Q is irreducible, then 

exp Q)t\ = (1 + o(l)) exp [t(7t, Q)t] wv T , 


elementwise as t oo, where w and v T are the left and right eigenvectors of A(ir,Q) corre¬ 
sponding to the dominant eigenvalue t(tt, Q) and normed so that v 7 w = 1. Since w and v are 
both strictly positive [ 3 ^, Theorem 2.6(b)], is the growth rate of the expected number 

of infected individuals during the initial stages of the epidemic. 

The following result shows that, under a certain condition on the recovery and infection rates, 
there is an optimal distribution of susceptible individuals which minimises the growth rate. 


Theorem 2.1. Let x{P,l) = (1 ~ Ejli TlMXEjli Pj X ) 1 and de fi ne K = (li + x{P, l))/Pi 
for i = 1,..., m. If 7 i > ~x(P, 7) for all i = 1,..., m, then 7 r* G <S and there exists a Q* G Q 
such that 

r(Ti-, Q*) > = t(tt*,Q) = x(P, 7), ( 2 - 4 ) 

for all it G S and all Q G Q. 


0 , 


Proof. The final equality in (12.411 holds as, for any Q G Q, A(yr*,Q)l = y(/3,7)l. From Seneta 
{ 31 ], Corollary 3, page 52], t(tt*,Q) = x(P,l)- 

To prove the inequality in (12.411 . fix Q G Q. As S is a convex set, it follows from Friedland 
Theorem 4.1] that t(-,Q) is a strictly convex functional on S . Therefore, 7 t(Q) will minimise 
r(-, Q) if and only if 

5r(vr,Q) 


E 

i=l 


dni 


7 r= 7 r (Q) 


(7Tj - 7T (Q)i) > 0, 


for all 7r G S. The partial derivatives of r(-,Q) are 

dT(ir,Q) fdiWiVi 


d-Ki 


w T v 


where w and v are the left and right eigenvectors of A(ir, Q ) corresponding to the dominant 


eigenvalue [ 31 , pg 183]. As [fn* — 7 \ = xi.P-,1 ) for i = 1,..., m, the eigenvectors of A(tt*,Q ) 
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and Q coincide so 1 is a right eigenvector of A{ n*, Q ) for any Q G Q. Therefore, the inequality 
in (12.41) will follow if there exists a Q* with left eigenvector w such that 

m 

Y, PiWi (77 - 7 r*) > 0, (2.5) 

i=l 

for all 7r G S. Set Q* such that QT = f 3 j l for i 7^ j. The left eigenvector of Q* satisfies Wi oc fif 1 . 
Therefore, inequality (|2.5D holds which proves t(tt,Q*) > t(tt*,Q*) for all 7 r € S. □ 

Under the condition of Theorem ETJ 7r* is the distribution of susceptible individuals which 
minimises the expected growth rate of the epidemic when infected individuals move to maximise 
the expected growth rate of the epidemic. A corresponding optimal migration rate matrix for 
susceptibles can easily be determined, but the optimal migration rate matrix is not unique. For 
example, any transition rate matrix R satisfying the detailed balance equations Rijn* = Rji^j, 
for all i,j G {1,. .. ,m} has 7 r* as its equilibrium distribution. Many other constructions are 
possible. 

Theorem 12.11 excludes the boundary case where 7 i > — x(/3, 7 ) for all i and 7 j = —x(/3,7) for 
some j. The difficulty in this case is that 7 r* 0 S since n* = 0 and there is no corresponding 
irreducible migration rate matrix for susceptibles. However, if we define, for any e > 0, 7 r e G S 
by 7 r| = (71 + e + x(/3, 7 )) / /3* for all i, then the calculations from the proof of Theorem 12.11 shows 
that (n e ,Q*) is an e-saddle point for 

The condition imposed in Theorem 12.11 will be satisfied provided the recovery rates do not 
vary too much between patches. In particular, if 7 * = 7 for i = , then x(/^, 7 ) = 

(Ylj= 1 ~ 7 > an d the condition holds for any /%, i = 1 ,... ,m. Also, if x(/5, 7 ) > 0 , then 

the condition must hold as the recovery rates are all nonnegative. On the other hand, when the 
condition does not hold, x(/?, 7 ) < 0. The following corollary shows that in this case there is a 
distribution of susceptible individuals for which the expected growth rate is negative, regardless 
of the movements of the infected individuals. 


Corollary 2.2. Suppose < 0. Then inf^s supg e Q t ( tt , Q) < 0. 

Proof. Define 77 = ^) _1 - Then for a11 Q G 2’ Q) 1 )* = 1 ~ !)> 1 = 

1 Now x(/ 3 j 7) < 0 implies ~^)~ l < 1- By Seneta [30, Corollary 3, page 52], 

t( if , Q) < 0 for all Q G Q. □ 


2.2. Optimising Rq, minor outbreak probability, and expected total size. When the 
condition of Theorem 12.11 does not hold, an alternative approach to controlling the disease 
spread is needed. As noted in Clancy % Section 4.1], the total size of the branching process 
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approximating the epidemic is the same as the total of an embedded Galton-Watson process. 
The behaviour of this Galton-Watson process is largely determined by the expected number of 
infectives produced by a single infective before its recovery. Denote by Aij(ir,Q) the expected 
number of infectives produced in group j by an individual first infected in group i. Then from 
Clancy j3, Section 4.1], A(7r, Q) = L(Q) diag(/3)diag(7r), where Lij{Q) is the expected amount of 
time that an individual who is first infected while in group i spends in group j before recovery. 

26, Proposition 2], L(Q) = (diag(7) — Q) -1 . Therefore, 


By Pollett and Stefanov 


A(vr, Q) = (diag(7) - Q) 1 diag(/3)diag(vr). 


The basic reproduction rate is the spectral radius of A(-7 t, Q ), which is denoted by Ro(tt, Q ). It 
is known that if Ro(tt, Q) < 1, then the Galton-Watson process goes extinct in finite time with 
probability one. Minimising Q ) provides an alternate means of limiting the spread of the 

disease. 

Theorem 2.3. Let w(/3,7) = (X^j=i T") 1 am ^ define iti = ^u;(/3,7) for i = 1 ,m. There 
exists a Q £ Q such that 


Ro(n, Q) > Rq (if, Q) = R 0 (fc, Q ) = u(/3, 7), 
for all tt € S and all Q € Q. 


( 2 . 6 ) 


Proof. The proof is similar to the proof of Theorem 12.11 We first prove the final equality in 
(12.61) . For any Q £ Q, Q1 = 0 so (diag(7) — Q)~ X 7 = 1. As 


A(tt,Q )1 = cj(/3,7) (diag(7) - Q) ^^(7)1 =cj(/3,7)1 


(2.7) 


it follows that Ro(n,Q) = cj(^, 7) 30|, Theorem 1.6]. 


To prove the inequality in (12.61) . fix Q G Q. As S is a convex set, Friedland [ill . Theorem 4.3] 
shows that Rq(-,Q) is a strictly convex functional on S. Therefore, tt(Q) minimises Rq(-,Q) if 
and only if 

dR 0 (TT,Q ) 


£ 

i =1 


dni 


T=ir(Q) 


(7 Ti ~ 7f(Q)j) > 0, 


for all 7r € S. The partial derivatives of Ro(-,Q ) are 

dRfi-n, Q) Pfi w (diag(q) - Q) -1 )*^ 


d-Ki 


w T v 


( 2 . 8 ) 


(2.9) 


where w and v are the left and right eigenvectors of A(n, Q ) corresponding to the dominant 
eigenvalue j^l, pg 183]. As noted previously, 1 is a right eigenvector of A (77, Q ) corresponding to 
the dominant eigenvalue uj(j3, 7). Substituting 7r = 7r and v = 1 in equation (12.91) and combining 









with equation (12.81) . we see that the inequality in (12.61) will follow if there exists a Q £ Q such 
that for all tt € S 

m 

^/3i(w(diag(7) - Q) _1 )i -*i) > 0. (2.10) 

i— 1 

As in Theorem EH set Q such that Qij = f3- 1 for i ^ j. Then the left eigenvector of A(fr,Q) 
corresponding to the dominant eigenvalue u(/3,'y) satisfies Wi oc 7 *//?* so (w(diag( 7 ) — oc 

/3~ ] . Therefore, inequality (|2.10l) holds which proves Rq(tt,Q) > Rq(tt,Q) for all tt G S. □ 


Theorem 12.31 shows that if is the distribution of susceptible individuals which minimises the 
basic reproduction rate of the epidemic when infected individuals move to maximise the basic 
reproduction rate of the epidemic. 

Although Theorems 12.11 and 12.31 give two different strategies for minimising the disease spread, 
the quantities x(/3, 7 ) and c o(/3, 7 ) are closely related. First, x(/3, 7 ) = 0 if and only if w(/3, 7 ) = 1 
in which case n* = if. Note we also have n* = tt when the recovery rates do not depend on 
the group. Second, < 0 if and only if uj(/3, 7 ) < 1. Therefore, Theorems 12.11 and 12.31 

and Corollary 12.21 imply that the i?o-optimal strategy yields R.q > 1 if and only if the r-optimal 
strategy yields r > 0 . 

Taking Rq as the objective function has the advantage that it is always possible to give an 
explicit optimal strategy. Although Rq is more tractable than r, it is only useful if the resulting 
optimal strategy reduces the extent of the original epidemic in some sense. In the context of 
ODE models, Diekmann et al [lo| notes that epidemics with high Rq do not necessarily have 
a fast increase of incidence. Therefore, one might question the relevance of reducing Rq if the 
threshold cannot be achieved and, if the threshold can be achieved, the advantage of reducing 
Rq further. To see why it is always useful to reduce Rq, it is necessary to consider the two cases 
cj(/ 3 , 7 ) < 1 and w(/3, 7 ) > 1 separately. The case where u(/3,'y) = 1 is not considered since in 
that case tt = tt* . 

We first examine how the optimal strategy from Theorem 12.31 relates to the probability 
that the branching process goes extinct in finite time. This probability is determined by the 
smallest fixed point of the probability generating function for the offspring distribution. For 
the branching process determined by (EH) EH, this probability generating function is 


<?i(u; 7T, Q) 


Qij u 3 + Pm u i + 7 i 
Qij + PiT^i T 7* 


The function g(u; tt, Q) = (gi(u; tt, Q),... , g m (u] tt, Q)) always has a fixed point at 1, that is 


g(l;7r,< 


= 1 . If t(tt,Q) > 0, then g(-;7r,Q) has a second fixed point, which is unique in 


(0, l) m Section 2.3]. Denote the smallest fixed point of g(-;7r, Q) in [0, l] m by q(7r, Q). The 
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probability of extinction in finite time is given by 

m 

lim P (Yj(t) = 0 for j = 1 ,..., m | 1 )( 0 ) = y* for i = 1 ,..., m) = TT qf{ n, Q). 

£—>•00 -*■ 

i= 1 

The following results shows that taking the distribution of susceptibles to be %r maximises the 
probability of extinction in finite time minimised over the starting location of the initial infected 
individual. 


Theorem 2.4. Ifu(/3, 7 ) > 1, then for any e > 0 there exists a Q e € Q such that 

mm{qi(Tr,Q e )} - e < min {qi(n, Q e )} = min {qi(n, Q)} = cj (/ 3 , 7) _1 (2.11) 

i i i 

for all 7 T S S and Q G Q. 


Proof. When it = 7r, the probability generating function of the offspring distribution is 

, _ +7< 

gi u;tt,Q = —- - —-- -. 

Ejy* Qij + 7 MA 7 ) +7i 

It can be verified by substitution that qi(n,Q) = u;(/3, 7 ) —1 , i = 1,..., m, for all Q G Q. It 
remains to prove the inequality in (12.1111 . which we achieve by determining an upper bound on 
q(7r, Q) for certain Q. 

As gi(-; 7 r, Q ) is a monotone function and q(7r, Q) is a fixed point of g(-; 7 r, Q), it follows that 
if, for some p € (0, l) m , g*(p; 7r ,Q) < Pi, i = 1,... , m, then qi(ir, Q) < Pi, i = 1,..., m. Let 
Si = YLj^i Qij an d choose Q € such that $i/Pi < e - For an y P e l) m > 


ffi(p;7r,Q e ) < 


Si + PmPj + 7 j 
Si + Pm + 7 i 


Si + + 7i ~ {Pm + 7i + 

+ PiTTi + 7i 


( 2 . 12 ) 


It can be verified by substitution into (12.1211 that 

'li + 


P* = 


is an upper bound on q( 7 r, Q e ). Therefore, 


Pm 


A 1 


min qi(TT,Q e ) — e < min 


7 i + & 

Pm 


— e. 


Suppose that, for some 7 r and all i = 1,..., m, 


7 i + Si 
Pm 


e > u{P, 7 ) \ 


then 


iTi < if i + 


SjUj(P, 7 ) 
Pi 


eTTiUliP,^). 


(2.13) 
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By summing over i in inequality (12.131) . we arrive at the contradiction 'Yffff 1 5i/ /3i > e. Therefore, 
for all 7 r € <S, there is at least one i such that 


H + 5j 

/3m 


which proves the inequality in ( 12 . 111 ) . 


e < w(/3,7) \ 


□ 


The previous theorem provides support for minimising Rq when cj(/3 ,7) > 1; it remains to 
justify minimising Rq when oj(/3, 7 ) < 1. Let Rj (tt, Q ) denote the number of individuals infected 
in node j starting from a single infected individual at node i. We now consider the effect of 
migration on the expected total size of the epidemic, 

r 7'( 7r > ( 9)j • 

The next result shows that taking the distribution of susceptibles to be tt minimises the total 
size of the epidemic maximised over the starting location of the initial infected individual. 


max IE 


Theorem 2.5. If u(/3, r )) < 1, then for any e > 0 there exists a Q e € Q such that 

/ m \ / rn 

maxE I Q e ) I + e > maxE I Q e 


(2.14) 


1.1=1 


^ j=i 


max E ( ^ TijifU, Q) ) = (1 — oj(/3, 7 )) 
7 =1 


-1 


for all tt € S and Q G Q. 


Proof. As the total size of the branching process approximating the epidemic is the same as that 
of an embedded Galton-Watson process whose offspring distribution has mean matrix A( 7 r, Q) 


a 


Section 4.1], 


E ( ] = (^A(7r,Q) ? 

G=> 


vr=0 


which is finite if and only if the spectral radius of A(n, Q ) is strictly less than one. From 
equation (12.71) . 

E =(1-o;(/3,7))- 1 , 

for all i = 1,..., m and all Q € Q. This proves the equality in (|2.14|) . To complete the proof, 
it remains to show that for any e > 0 there exists a Q e € Q, such that 


maxE j Rj (tt , Q e ) ) +e > (1 -w(^,7)) \ 
U =1 


(2.15) 
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for all 7r £ S. Take Q e = d/ll 1 — ml) where 5 satisfies 

m 

0 < dm^2 7;" 1 < e (l - w(/3,7))(o;(/3,7) _1 - 1). 

i —1 

Applying the Woodbury matrix identity to A(ir, Q e ), we obtain 


A(vr ,Q £ ) 



5 

(l-dl^Fj 1 !) 



r< 5 1 diag(/3)diag(vr), 


(2.16) 


where T^ = diag(y + 5ml). Hence, in the partial order of positive matrices, 


diag(/3)diag(7r)T (5 1 < A(n,Q e ). 


Therefore, the expected total size is finite for all e > 0 only if <1, i = 1,... ,m, in 

which case 

— 1 / m 


l - 


Pm 


7 i + 5m 


<E Y, T *M) » 


(2.17) 


J =i 


for alii = 1,..., m. Suppose that inequality ([2.151) did not hold for some it £ S. Then, for all 

h 


E ( ) +e < (i-w(^,7)) 1 . 

U =1 


Inequality (12.171) would then implies 


1 - 


Pm 


7 i + 5m 

This inequality can be rearranged to 


TTi - e(l -w(/3,7))^ ( 1 - ^ 


+ e < (1 - w(^, 7 )) 


-i 


Pi \ 7i + 


. 5m 

> ( 1 -I TU. 

7 j + ora 


Summing this inequality over i, we find 


e(l- 


m 


UJ 


t=l 


a \ m 

) <gm y 
7 ,- + dm ' 


TTi 


As 5 is chosen to satisfy the inequality (12.161) . we obtain a contradiction. Hence, inequality 
(12.151) holds for all it £ S. □ 
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3. Numerical comparisons 


In this section we investigate numerically two issues. The first issue concerns the optimal 
distribution of susceptibles with respect to minimising the expected growth rate. Theorem 12.11 
gives the optimal distribution only if 7, > — f° r all *■ We have also seen that if 7j > 

— x(( 3 , 7 ) for all i and 7 j = — x(/ 3 , 7 ) for some j, then there is a sequence of distributions Tr e which 
achieve within e the optimal value of r and for which lim e _>.o vr| = 0. We might expect that if 
the recovery rate for this group were to decrease, then the optimal distribution would place no 
susceptibles in group j. This was investigated in a two group epidemic with fi\ = 1, fa = 2, and 
71 and 72 in (0,4). For these epidemics both supggg infogs t ( tt , Q) and infogs supggo t ( it , Q) 
were computed by nested optimisation using the optim and optimize functions in R [27]. The 


two quantities differed by less than 10~ 4 in all instances computed. The optimal value of r(7r, Q) 
is plotted in Figure [0 Note that in most of the region plotted growth rate is negative. This is 
to be expected as when the condition of Theorem 12.11 does not hold, the optimal growth rate 
must be negative from Corollary 12.21 The numerical results confirms our intuition that tt* = 0 
if 7 i < — x(P>l) fo r the two group model. From the plot it is seen that if 71 < — x(/3,7) < 72, 
then increasing 72 has no effect on the optimal value of t(tt,Q). This is explained as when 
71 < —x(/3,7) < 72, the optimal distribution of susceptibles has tt \ = 1. On the other hand, 
the inequality 71 < — x(/3,7) implies 71 < 72 — ^2, so the infectives slow the decrease of the 
epidemic by moving to group one. Therefore, increasing the recovery rate in group two has no 
effect on the growth rate of the epidemic and infogs supg g g t(tt, Q) = —71. 

By construction, tt* and if are the optimal distribution of susceptibles for minimising r and Rq 
respectively. We now consider their performance on the alternate criteria, that is we calculate 
supggQ Rq(tt*, Q) and supg g g r(7r, Q). Figures [2] and [3] show how much these quantities are 
increased by taking alternate optimal distributions of susceptibles. Qualitatively, both figures 
are very similar. Both quantities plotted achieve their minimum for the same set of 7, indicated 
by the dashed line, as tt* = n for these values of 7. Also, an abrupt change in the contours 
occur along the dashed line in both figures. This is due to n* placing zero probability in one of 
the groups for those values of 7 outside the dashed lines. 

For most values of 7 the optimal choice for one criterion appears to result in reasonable 
performance in the other. In particular, in the region where 7! « 72 « + P 2 ), t 

and Rq take approximately the same value under -k* and tt. However, for small values of 7, 
the performance of the alternate distributions rapidly deteriorates for both r and Rq. This is 
expected as when 7 is small, x(/3,7) tends to be large which causes the difference between tt* 
and tt to also be large. 
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optimal t(ji, Q) 



Yi 


Figure 1. A contour plot of the optimal value of for the two group 

model with f 3 \ = 1 and @2 = 2. The region between the dotted lines corresponds 
to the region where 7j > —x(/3,7) for i = 1,2. Note that where the contour lines 
are horizontal, t(tt, Q ) = —72 and where the contour lines are vertical t(tt, Q ) = 
-7i- 


4. Discussion 

The conclusions of Theorems m and 12.31 are in part not surprising; in order to minimise 
the spread of the disease most susceptible individuals should belong to groups with relatively 
low infection rates and high recovery rates. However, for the form of contact rate assumed 
here, this needs to be balanced with the fact that contact rates are higher in groups with 
larger populations. Although Theorems 12.41 and 12.51 showed that this balance is achieved in 
the same way for i?Q-optimal, extinction probability optimal, and expected total size optimal 
distributions of susceptibles, it was achieved differently for r-optimal distribution of susceptibles. 
It is conceivable that this balance might be achieved differently for other measures of disease 
spread. 

In our analysis, we have focussed on the branching process approximation to the epidemic. 
Another widely used approximation is provided by the solution to an ordinary differential equa¬ 
tion (ODE). Assume that infected individuals recover without immunity. For the epidemic 
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log(sup R 0 (ji*. Q)/to(P, y)) 



Yi 


Figure 2. A contour plot of the optimal value of log(supQ g g Rq(7t* , Q)/u(f3, 7 )) 
for the two group model with B\ = 1 and ,02 = 2. The dashed lines correspond to 
where u(/3,'y) = 1 so n* = n. The region between the dotted lines corresponds 
to the region where 7 j > — x(/3,7) for i = 1 , 2 . 


described at the beginning of Section [2l the ODE approximation is given by the solution to 


dxi 

dt 


T, RjiXjit) + RaXiit) + 7 iyi(t) - /3iXi(t)yi(t ) 
j¥=i 


dyi 

dt 


y QjiVj{t) + Qayiit) - 7 iVi(t) + PiXi(t)yi(t). 


(4.18) 

(4.19) 


It is known that if N 1 2Q(0) A 27 ( 0 ) and N 1 T)(0) A y*(0) for i = 1,..., m as N —> 00, then 
for any finite T > 0 and any e > 0 


lim Pr ( sup ( Y] \N 1 X i (t) - Xi(t)\ + V \N l Yi{t) - yi(t)\ ) >e) = 0, 

N -+°° \te[°,T] ^ ' 1 


i =1 


see 


21 


Theorems 12.11 and 12.31 can still be used to determine the optimal distribution of susceptibles 
for the ODE model (14.1811 - (14.191) . First, consider the application of Theorem 12.11 The spectrum 
of the Jacobian of the ODE model at the disease free equilibrium is given by the union of the 
spectrum of A(tt,Q), where 7 r is the unique solution to ttR = 0 subject to 7 rl = 1, and the 
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supt(S, Q)-supt(ji‘, Q) 



Yi 

Figure 3. A contour plot of the optimal value of supg g g t(tt, Q) — 
supg g Q t( 7 t*, Q) for the two group model with /?i = 1 and /?2 = 2. The dashed 
lines correspond to where x(Pil) = 0 so it* = if. The region between the dotted 
lines corresponds to the region where 7 j > —x(/3, 7 ) for i = 1 ,2. 

spectrum of R with the zero eigenvalue removed. Therefore, if 7 j > —x{Pil) for i = 1, ■ ■ ■ ,m, 
then Theorem 12.II determines the r optimal choice of n. However, for this to be attained, R must 
be chosen so that the non-zero eigenvalues of R have real part less than x(0, 7 ). Theorem 12.31 
can similarly be applied to the ODE model. The next generation matrix [331 . Section 3] for 
the ODE model is given by A(tt,Q) t . As the basic reproduction number for the ODE model 
is given by the spectral radius of the next generation matrix, Theorem 12.31 determines the Rq 
optimal distribution of susceptibles in the metapopulation. We are unaware of an interpretation 
of Theorems 12.41 and 12.51 for the ODE model. 

We have previously noted that the desire for susceptible individuals to belong to a group 
with a low infection rate and high recovery rate needs to be balanced with the fact that contact 
rates are higher in groups with larger populations. This was due to the assumption of density 
dependent contact rates. An alternative is to assume frequency dependent contact rates, that 
is to assume the per capita contact rate in a group does not depend on the size of the group. 
Allen et al { 2 ] studied a frequency-dependent SIS metapopulation model, which in our notation 
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is given by 


dxi 

dt 

dyi 

dt 


y, RjiXjit) + RaXi(t) + 7 iViit) - 
i¥=i 

y QjiVj{t) + QnViit) - ny^t) + 

i¥=i 


/%Xj(t)yj(t ) 
Xi(t) +yi{t) 

PiXi(t)yj(t) 
Xi{t ) + yi(t) ' 


For this model, the next generation matrix is given by diag(/3)(diag(y) — Q T ) -1 [2|, Lemma 
2.2] so Rq does not depend on the migration rates of susceptible individuals. Therefore, we 
are unable to control the disease spread through the altering the migration rates of susceptible 
individuals. Although frequency dependent and density dependent contact rates are the most 
commonly assumed form for contact rates, it is possible to consider contact rates that are some 
general function of the size of the group. For these more general contact rates, we expect results 
similar to Theorems 12.11 and 12.31 to hold. 
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